#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>

static int _is_same(const cv::Mat& m1, const cv::Mat& m2) {
    cv::Mat diff = m1 != m2;
    std::vector<cv::Point> nonzeros;
    cv::findNonZero(diff, nonzeros);
    return nonzeros.size();
}

static void _resize_nearest(const cv::Mat& src, cv::Mat& dst, double fx, double fy) {
    for (int i = 0; i < dst.rows; i++) {
        uchar* pd = dst.ptr<uchar>(i);
        int sy = std::min(cvFloor(i/fy), src.rows-1);
        const uchar* ps = src.ptr<uchar>(sy);
        for (int j = 0; j < dst.cols; j++) {
            int sx = std::min(cvFloor(j/fx), src.cols-1);
            pd[j] = ps[sx];
        }
    }
}

static int _clip(int x, int a, int b) {
    if (x >= a) {
        if (x < b) {
            return x;
        }
        return b-1;
    }
    return a;
}

static void _resize_bilinear(const cv::Mat& src, cv::Mat& dst, double x_dst_2_src, double y_dst_2_src, bool isArea=false) {
    std::vector<int> x_of_src(dst.cols);
    std::vector<int> y_of_src(dst.rows);
    std::vector<short> alpha0(dst.cols), alpha1(dst.cols);
    std::vector<short> beta0(dst.rows), beta1(dst.rows);
    
    int ksize = 2;
    int ksize2 = ksize / 2;
    int xmin = 0, xmax = dst.cols;
    
    for (int dx = 0; dx < dst.cols; dx++) {
        float fx;
        int sx;
        if (!isArea) {
            fx = (float)((dx+0.5) / x_dst_2_src - 0.5);
            sx = cvFloor(fx);
            fx -= sx;
        } else {
            sx = cvFloor(dx / x_dst_2_src);
            fx = (float)((dx+1)-(sx+1)*x_dst_2_src);
            fx = (fx <= 0 ? 0.f : fx-cvFloor(fx));
        }
        if (sx < 0) {
            xmin = dx+1;
            fx = 0;
            sx = 0;
        }
        if (sx >= src.cols-1) {
            xmax = std::min(xmax, dx);
            fx = 0;
            sx = src.cols-1;
        }
        x_of_src[dx] = sx;
        float u0 = 1.f - fx;
        float u1 = fx;
        alpha0[dx] = cv::saturate_cast<short>(u0 * 2048);
        alpha1[dx] = cv::saturate_cast<short>(u1 * 2048);
    }
    
    for (int dy = 0; dy < dst.rows; dy++) {
        float fy;
        int sy;
        if (!isArea) {
            fy = (float)((dy+0.5) / y_dst_2_src - 0.5);
            sy = cvFloor(fy);
            fy -= sy;
        } else {
            sy = cvFloor(dy / y_dst_2_src);
            fy = (float)((dy+1)-(sy+1)*y_dst_2_src);
            fy = (fy <= 0 ? 0.f : fy-cvFloor(fy));
        }
        y_of_src[dy] = sy;
        float v0 = 1.f - fy;
        float v1 = fy;
        beta0[dy] = cv::saturate_cast<short>(v0 * 2048);
        beta1[dy] = cv::saturate_cast<short>(v1 * 2048);
    }
    
    for (int dy = 0; dy < dst.rows; dy++) {
        int sy = y_of_src[dy];
        int sy0 = _clip(sy, 0, src.rows);
        int sy1 = _clip(sy+1, 0, src.rows);
        std::vector<int> D0(dst.cols), D1(dst.cols);
        const uchar* S0 = src.ptr<uchar>(sy0);
        const uchar* S1 = src.ptr<uchar>(sy1);
        int dx = 0;
        for (; dx < xmax; dx++) {
            int sx = x_of_src[dx];
            int a0 = alpha0[dx], a1 = alpha1[dx];
            int t0 = S0[sx]*a0 + S0[sx+1]*a1;
            int t1 = S1[sx]*a0 + S1[sx+1]*a1;
            D0[dx] = t0;
            D1[dx] = t1;
        }
        for (; dx < dst.cols; dx++) {
            int sx = x_of_src[dx];
            D0[dx] = int(S0[sx]*2048);
            D1[dx] = int(S1[sx]*2048);
        }
        
        short b0 = beta0[dy];
        short b1 = beta1[dy];
        for (int x = 0; x < dst.cols; x++) {
            //dst.ptr<uchar>(dy)[x] = uchar((b0*D0[x]+b1*D1[x]) >> 22);
            dst.ptr<uchar>(dy)[x] = uchar((((b0*(D0[x]>>4)) >> 16) + ((b1*(D1[x]>>4)) >> 16) + 2) >> 2);
        }
    }
}

static void _cubic_interpolate_coeffs(float x, float& c1, float& c2, float& c3, float& c4) {
    const float A = -0.75f;
    /*c1 = A*x*x*x - 2*A*x*x + A*x;
    c2 = (A+2)*x*x*x - (A+3)*x*x + 1;
    c3 = -(A+2)*x*x*x + (2*A-3)*x*x - A*x;
    //c4 = -A*x*x*x + A*x*x;
    c4 = 1.f - c1 - c2 - c3;
    */
    c1 = ((A*(x+1)-5*A) * (x+1) + 8*A) * (x+1) - 4*A;
    c2 = ((A+2)*x - (A+3))*x*x+1;
    c3 = ((A+2)*(1-x)-(A+3))*(1-x)*(1-x)+1;
    c4 = 1.f - c1 - c2 - c3;
}

static void _resize_cubic(const cv::Mat& src, cv::Mat& dst, double x_dst_2_src, double y_dst_2_src) {
    std::vector<int> x_of_src(dst.cols);
    std::vector<int> y_of_src(dst.rows);
    std::vector<short> alpha0(dst.cols), alpha1(dst.cols), alpha2(dst.cols), alpha3(dst.cols);
    std::vector<short> beta0(dst.rows), beta1(dst.rows), beta2(dst.rows), beta3(dst.rows);
    
    int ksize = 4;
    int ksize2 = ksize / 2;
    int xmin = 0, xmax = dst.cols;
    
    for (int dx = 0; dx < dst.cols; dx++) {
        float fx = (float)((dx+0.5) / x_dst_2_src - 0.5);
        int sx = cvFloor(fx);
        fx -= sx;
        if (sx < ksize2 - 1) {
            xmin = dx+1;
        }
        if (sx >= src.cols-ksize2) {
            xmax = std::min(xmax, dx);
        }
        x_of_src[dx] = sx;
        float c1, c2, c3, c4;
        _cubic_interpolate_coeffs(fx, c1, c2, c3, c4);
        alpha0[dx] = cv::saturate_cast<short>(c1 * 2048);
        alpha1[dx] = cv::saturate_cast<short>(c2 * 2048);
        alpha2[dx] = cv::saturate_cast<short>(c3 * 2048);
        alpha3[dx] = cv::saturate_cast<short>(c4 * 2048);
    }
    
    for (int dy = 0; dy < dst.rows; dy++) {
        float fy = (float)((dy+0.5) / y_dst_2_src - 0.5);
        int sy = cvFloor(fy);
        fy -= sy;
        y_of_src[dy] = sy;
        float c1, c2, c3, c4;
        _cubic_interpolate_coeffs(fy, c1, c2, c3, c4);
        beta0[dy] = cv::saturate_cast<short>(c1 * 2048);
        beta1[dy] = cv::saturate_cast<short>(c2 * 2048);
        beta2[dy] = cv::saturate_cast<short>(c3 * 2048);
        beta3[dy] = cv::saturate_cast<short>(c4 * 2048);
    }
    
    for (int dy = 0; dy < dst.rows; dy++) {
        int sy = y_of_src[dy];
        int sy0 = _clip(sy-1, 0, src.rows);
        int sy1 = _clip(sy, 0, src.rows);
        int sy2 = _clip(sy+1, 0, src.rows);
        int sy3 = _clip(sy+2, 0, src.rows);
        std::vector<int> D0(dst.cols), D1(dst.cols), D2(dst.cols), D3(dst.cols);
        const uchar* S0 = src.ptr<uchar>(sy0);
        const uchar* S1 = src.ptr<uchar>(sy1);
        const uchar* S2 = src.ptr<uchar>(sy2);
        const uchar* S3 = src.ptr<uchar>(sy3);
        int dx = 0;
        for (; dx < xmin; dx++) {
            int sx = x_of_src[dx];
            int sx0 = _clip(sx-1, 0, src.cols);
            int sx1 = _clip(sx, 0, src.cols);
            int sx2 = _clip(sx+1, 0, src.cols);
            int sx3 = _clip(sx+2, 0, src.cols);
            D0[dx] = S0[sx0]*alpha0[dx]+S0[sx1]*alpha1[dx]+S0[sx2]*alpha2[dx]+S0[sx3]*alpha3[dx];
            D1[dx] = S1[sx0]*alpha0[dx]+S1[sx1]*alpha1[dx]+S1[sx2]*alpha2[dx]+S1[sx3]*alpha3[dx];
            D2[dx] = S2[sx0]*alpha0[dx]+S2[sx1]*alpha1[dx]+S2[sx2]*alpha2[dx]+S2[sx3]*alpha3[dx];
            D3[dx] = S3[sx0]*alpha0[dx]+S3[sx1]*alpha1[dx]+S3[sx2]*alpha2[dx]+S3[sx3]*alpha3[dx];
        }
        for (; dx < xmax; dx++) {
            int sx = x_of_src[dx];
            D0[dx] = S0[sx-1]*alpha0[dx]+S0[sx]*alpha1[dx]+S0[sx+1]*alpha2[dx]+S0[sx+2]*alpha3[dx];
            D1[dx] = S1[sx-1]*alpha0[dx]+S1[sx]*alpha1[dx]+S1[sx+1]*alpha2[dx]+S1[sx+2]*alpha3[dx];
            D2[dx] = S2[sx-1]*alpha0[dx]+S2[sx]*alpha1[dx]+S2[sx+1]*alpha2[dx]+S2[sx+2]*alpha3[dx];
            D3[dx] = S3[sx-1]*alpha0[dx]+S3[sx]*alpha1[dx]+S3[sx+1]*alpha2[dx]+S3[sx+2]*alpha3[dx];
        }
        for (; dx < dst.cols; dx++) {
            int sx = x_of_src[dx];
            int sx0 = _clip(sx-1, 0, src.cols);
            int sx1 = _clip(sx, 0, src.cols);
            int sx2 = _clip(sx+1, 0, src.cols);
            int sx3 = _clip(sx+2, 0, src.cols);
            D0[dx] = S0[sx0]*alpha0[dx]+S0[sx1]*alpha1[dx]+S0[sx2]*alpha2[dx]+S0[sx3]*alpha3[dx];
            D1[dx] = S1[sx0]*alpha0[dx]+S1[sx1]*alpha1[dx]+S1[sx2]*alpha2[dx]+S1[sx3]*alpha3[dx];
            D2[dx] = S2[sx0]*alpha0[dx]+S2[sx1]*alpha1[dx]+S2[sx2]*alpha2[dx]+S2[sx3]*alpha3[dx];
            D3[dx] = S3[sx0]*alpha0[dx]+S3[sx1]*alpha1[dx]+S3[sx2]*alpha2[dx]+S3[sx3]*alpha3[dx];
        }
        
        
        short b0 = beta0[dy];
        short b1 = beta1[dy];
        short b2 = beta2[dy];
        short b3 = beta3[dy];
        for (int x = 0; x < dst.cols; x++) {
            dst.ptr<uchar>(dy)[x] = cv::saturate_cast<uchar>(((D0[x]*b0+D1[x]*b1+D2[x]*b2+D3[x]*b3) + (1 << 21)) >> 22);
        }
    }
}

static void _lanczos4_interpolate_coeffs(float x, float* coeffs) {
    /*if (x < FLT_EPSILON) {
        for (int i = 0; i < 8; i++) {
            coeffs[i] = 0;
        }
        coeffs[3] = 1;
        return;
    }
    float sum = 0;
    for (int i = 0; i < 8; i++) {
        double xi = -(x+3-i)*CV_PI;
        coeffs[i] = 4*std::sin(xi/4)*std::sin(xi)/(xi*xi);
        sum += coeffs[i];
    }
    for (int i = 0; i < 8; i++) {
        coeffs[i] /= sum;
    }*/
    const double s45 = 0.70710678118654752440084436210485;
    const double cs[][2] = {
        {1, 0}, {-s45, -s45}, {0, 1}, {s45, -s45}, {-1, 0}, {s45, s45}, {0, -1}, {-s45, s45}  
    };
    
    float sum = 0;
    double y0 = -(x+3)*CV_PI*0.25, s0 = std::sin(y0), c0 = std::cos(y0);
    for (int i = 0; i < 8; i++) {
        float y0_ = (x+3-i);
        if (fabs(y0_) >= 1e-6f) {
            double y = -y0_*CV_PI*0.25;
            coeffs[i] = (float)((cs[i][0]*s0+cs[i][1]*c0)/(y*y));
        } else {
            coeffs[i] = 1e30f;
        }
        sum += coeffs[i];
    }
    sum = 1.f / sum;
    for (int i = 0; i < 8; i++) {
        coeffs[i] *= sum;
    }
}

static void _resize_lanczos4(const cv::Mat& src, cv::Mat& dst, double x_dst_2_src, double y_dst_2_src) {
    std::vector<int> x_of_src(dst.cols);
    std::vector<int> y_of_src(dst.rows);
    int ksize = 8, ksize2 = ksize / 2;
    
    std::vector<std::vector<short>> alpha(ksize), beta(ksize);
    for (int k = 0; k < ksize; k++) {
        alpha[k].resize(dst.cols);
        beta[k].resize(dst.rows);
    }
    
    int xmin = 0, xmax = dst.cols;
    
    for (int dx = 0; dx < dst.cols; dx++) {
        float fx = (float)((dx+0.5) / x_dst_2_src - 0.5);
        int sx = cvFloor(fx);
        fx -= sx;
        if (sx < ksize2 - 1) {
            xmin = dx+1;
        }
        if (sx >= src.cols-ksize2) {
            xmax = std::min(xmax, dx);
        }
        x_of_src[dx] = sx;
        float cs[8];
        _lanczos4_interpolate_coeffs(fx, cs);
        for (int k = 0; k < ksize; k++) {
            alpha[k][dx] = cv::saturate_cast<short>(cs[k]*2048);
        }
    }
    
    for (int dy = 0; dy < dst.rows; dy++) {
        float fy = (float)((dy+0.5) / y_dst_2_src - 0.5);
        int sy = cvFloor(fy);
        fy -= sy;
        y_of_src[dy] = sy;
        float cs[8];
        _lanczos4_interpolate_coeffs(fy, cs);
        for (int k = 0; k < ksize; k++) {
            beta[k][dy] = cv::saturate_cast<short>(cs[k]*2048);
        }
    }
    
    std::vector<std::vector<int>> DRows(ksize);
    std::vector<int> sys(ksize);
    for (int k = 0; k < ksize; k++) {
        DRows[k].resize(dst.cols);
    }
    
    for (int dy = 0; dy < dst.rows; dy++) {
        int sy = y_of_src[dy];
        for (int k = 0; k < ksize; k++) {
            sys[k] = _clip(sy-ksize2+1+k, 0, src.rows);
        }
        
        for (int k = 0; k < ksize; k++) {
            const uchar* S = src.ptr<uchar>(sys[k]);
            std::vector<int>& D = DRows[k];
            int dx = 0;
            for (; dx < xmin; dx++) {
                int sx = x_of_src[dx]-3;
                D[dx] = 0;
                for (int j = 0; j < ksize; j++) {
                    int sxj = _clip(sx+j, 0, src.cols);
                    D[dx] += S[sxj]*alpha[j][dx];
                }
            }
            for (; dx < xmax; dx++) {
                int sx = x_of_src[dx];
                D[dx] = 0;
                for (int j = 0; j < ksize; j++) {
                    D[dx] += S[sx-3+j]*alpha[j][dx];
                }
            }
            for (; dx < dst.cols; dx++) {
                int sx = x_of_src[dx]-3;
                D[dx] = 0;
                for (int j = 0; j < ksize; j++) {
                    int sxj = _clip(sx+j, 0, src.cols);
                    D[dx] += S[sxj]*alpha[j][dx];
                }
            }
        }        
        
        short b[8];
        for (int ii = 0; ii < ksize; ii++) {
            b[ii] = beta[ii][dy];
        }
        for (int x = 0; x < dst.cols; x++) {
            int res = 0;
            for (int ii = 0; ii < ksize; ii++) {
                res += DRows[ii][x]*b[ii];
            }
            dst.ptr<uchar>(dy)[x] = cv::saturate_cast<uchar>((res + (1 << 21)) >> 22);
        }
    }
}

static int _compute_resize_area_coeffs(int ssize, int dsize, double dst_2_src, std::vector<int>& src_indices, std::vector<int>& dst_indices, std::vector<float>& alpha) {
    double scale = 1. / dst_2_src;
    for (int i = 0; i < dsize; i++) {
        double fs1 = i*scale;
        double fs2 = fs1 + scale;
        double cell_width = std::min(scale, ssize-fs1);
        int s1 = cvCeil(fs1);
        int s2 = cvFloor(fs2);
        s2 = std::min(s2, ssize-1);
        s1 = std::min(s1, s2);
        if (s1 - fs1 > 1e-3) {
            dst_indices.push_back(i);
            src_indices.push_back(s1-1);
            alpha.push_back((float)((s1-fs1)/cell_width));
        }
        for (int s = s1; s < s2; s++) {
            dst_indices.push_back(i);
            src_indices.push_back(s);
            alpha.push_back((float)(1.0/cell_width));
        }
        if (fs2 - s2 > 1e-3) {
            dst_indices.push_back(i);
            src_indices.push_back(s2);
            alpha.push_back((float)(std::min(std::min(fs2-s2, 1.0), cell_width) / cell_width));
        }
    }
    return alpha.size();
}

static void _resize_area(const cv::Mat& src, cv::Mat& dst, double x_dst_2_src, double y_dst_2_src) {
    double scale_x = 1. / x_dst_2_src, scale_y = 1. / y_dst_2_src;
    
    if (scale_x >= 1 && scale_y >= 1) {
        std::vector<int> src_x_indices, dst_x_indices, src_y_indices, dst_y_indices;
        std::vector<float> x_alpha, y_beta;
        int x_indices_size = _compute_resize_area_coeffs(src.cols, dst.cols, x_dst_2_src, src_x_indices, dst_x_indices, x_alpha);
        int y_indices_size = _compute_resize_area_coeffs(src.rows, dst.rows, y_dst_2_src, src_y_indices, dst_y_indices, y_beta);
        
        std::vector<float> sum(dst.cols), buf(dst.cols);
        for (int i = 0; i < dst.cols; i++) {
            sum[i] = 0;
        }
        for (int dy = 0; dy < dst.rows; dy++) {
            for (int j = 0; j < y_indices_size; j++) {
                if (dy == dst_y_indices[j]) {
                    int sy = src_y_indices[j];
                    float beta = y_beta[j];
                    
                    const uchar* S = src.ptr<uchar>(sy);
                    for (int dx = 0; dx < dst.cols; dx++) {
                        buf[dx] = 0;
                    }
                    for (int i = 0; i < x_indices_size; i++) {
                        int dx = dst_x_indices[i];
                        int sx = src_x_indices[i];
                        float alpha = x_alpha[i];
                        buf[dx] += S[sx]*alpha;
                    }
                    for (int dx = 0; dx < dst.cols; dx++) {
                        sum[dx] += beta*buf[dx];
                    }
                }
            }
            for (int dx = 0; dx < dst.cols; dx++) {
                dst.ptr<uchar>(dy)[dx] = cv::saturate_cast<uchar>(sum[dx]);
                sum[dx] = 0;
            }
        }
    } else {
        _resize_bilinear(src, dst, x_dst_2_src, y_dst_2_src, true);
    }
}

static cv::Mat _my_resize(const cv::Mat& src, cv::Size dsize, double fx=0, double fy=0, int interpolation=cv::INTER_LINEAR) {
    CV_Assert(src.type() == CV_8UC1);
    
    cv::Size ssize = src.size();
    CV_Assert(!ssize.empty());
    
    if (dsize.empty()) {
        CV_Assert(fx > 0 && fy > 0);
        dsize = cv::Size(cv::saturate_cast<int>(ssize.width*fx), cv::saturate_cast<int>(ssize.height*fy));
        CV_Assert(!dsize.empty());
    } else {
        fx = (double)dsize.width / ssize.width;
        fy = (double)dsize.height / ssize.height;
        CV_Assert(fx > 0 && fy > 0);
    }
    
    cv::Mat dst(dsize, src.type());
    if (dsize == ssize) {
        src.copyTo(dst);
        return dst;
    }
    
    if (interpolation == cv::INTER_NEAREST) {
        _resize_nearest(src, dst, fx, fy);
        return dst;
    }
    
    if (interpolation == cv::INTER_LINEAR) {
        _resize_bilinear(src, dst, fx, fy);
        return dst;
    }
    
    if (interpolation == cv::INTER_CUBIC) {
        _resize_cubic(src, dst, fx, fy);
        return dst;
    }
    
    if (interpolation == cv::INTER_LANCZOS4) {
        _resize_lanczos4(src, dst, fx, fy);
        return dst;
    }
    
    if (interpolation == cv::INTER_AREA) {
        _resize_area(src, dst, fx, fy);
        return dst;
    }
    
    return dst;
}

int main(int argc, char **argv) {
    cv::Mat src = cv::imread("/home/xlll/Downloads/opencv/samples/data/lena.jpg", cv::IMREAD_GRAYSCALE);
    if (src.empty()) {
        std::cout << "failed to read image" << std::endl;
        return EXIT_FAILURE;
    }
    
    cv::Size dsize1(src.cols / 2, src.rows * 2 / 3);
    cv::Size dsize2((int)(src.cols*1.2), (int)(src.rows*1.5));
    cv::Mat dst1, dst2;
    cv::Mat dst11, dst22;
    //cv::resize(src, dst1, dsize1);
    //cv::resize(src, dst2, dsize2);
    
    /*
    cv::resize(src, dst1, dsize1, 0, 0, cv::INTER_NEAREST);
    cv::resize(src, dst2, dsize2, 0, 0, cv::INTER_NEAREST);
    
    dst11 = _my_resize(src, dsize1, 0, 0, cv::INTER_NEAREST);
    dst22 = _my_resize(src, dsize2, 0, 0, cv::INTER_NEAREST);
    */
    
    /*
    cv::resize(src, dst1, dsize1, 0, 0, cv::INTER_LINEAR);
    cv::resize(src, dst2, dsize2, 0, 0, cv::INTER_LINEAR);
    
    dst11 = _my_resize(src, dsize1, 0, 0, cv::INTER_LINEAR);
    dst22 = _my_resize(src, dsize2, 0, 0, cv::INTER_LINEAR);
    */
    /*
    cv::resize(src, dst1, dsize1, 0, 0, cv::INTER_CUBIC);
    cv::resize(src, dst2, dsize2, 0, 0, cv::INTER_CUBIC);
    
    dst11 = _my_resize(src, dsize1, 0, 0, cv::INTER_CUBIC);
    dst22 = _my_resize(src, dsize2, 0, 0, cv::INTER_CUBIC);
    */
    /*
    cv::resize(src, dst1, dsize1, 0, 0, cv::INTER_LANCZOS4);
    cv::resize(src, dst2, dsize2, 0, 0, cv::INTER_LANCZOS4);
    
    dst11 = _my_resize(src, dsize1, 0, 0, cv::INTER_LANCZOS4);
    dst22 = _my_resize(src, dsize2, 0, 0, cv::INTER_LANCZOS4);
    */
    /*
    cv::resize(src, dst1, dsize1, 0, 0, cv::INTER_AREA);
    cv::resize(src, dst2, dsize2, 0, 0, cv::INTER_AREA);
    
    dst11 = _my_resize(src, dsize1, 0, 0, cv::INTER_AREA);
    dst22 = _my_resize(src, dsize2, 0, 0, cv::INTER_AREA);
    */
    
    cv::Size dsize(src.cols/3, src.rows/2);
    cv::resize(src, dst1, dsize, 0, 0, cv::INTER_LINEAR);
    cv::resize(src, dst11, dsize, 0, 0, cv::INTER_AREA);
    
    cv::resize(src, dst2, dsize, 0, 0, cv::INTER_LINEAR_EXACT);
    cv::resize(src, dst22, dsize, 0, 0, cv::INTER_AREA);
    
    int nonzerosize1 = _is_same(dst1, dst11);
    int nonzerosize2 = _is_same(dst2, dst22);
    cv::String str1 = (nonzerosize1 == 0 ? "Same" : "Not Same: " + std::to_string(nonzerosize1));
    cv::String str2 = (nonzerosize2 == 0 ? "Same" : "Not Same: " + std::to_string(nonzerosize2));
    std::cout << str1 << std::endl;
    std::cout << str2 << std::endl;
    
    cv::imshow("src", src);
    cv::imshow("dst1", dst1);
    cv::imshow("dst2", dst2);
    cv::waitKey(0);
    return EXIT_SUCCESS;
}
